"Kari Kivisalo" <kar### [at] kivisalonet> wrote in message
news:3995801A.AC206C81@kivisalo.net...
> Chris Huff wrote:
> >
> > Could you give more details on how you calculate
> > it? I didn't see an in-depth explanation on the site.
>
> Here is the theory on 2d electrical potential and how to obtain numerical
> approximation using difference equations and iterations. This can be
solved
> with FFT also but I don't yet know how. Should be easy enough to find:
> "solving poisson with fft". Numerical Recipes has nice code for 3D FFT
> which calculates the solution in-place saving memory.
>
> Laplace: u(xx) + u(yy) = 0
> u(xx) = 2nd derivative of u over variable x
>
> Here u is potential (hf height). x and y are hf coordinates.
> Resulting difference formula in free areas (x,y) + 1 pixel boundary:
>
> u(x+1,y) + u(x,y+1) + u(x-1,y) + u(x,y-1) = 4*u(x,y)
>
> So basically the procedure is: Take surrounding 4 pixels and put
> the average in the middle. Repeat for all free pixels (x,y). If
> there are n free pixels the resulting linear system is n*5. I solve
> this with Gauss-Seidel iteration.
>
> More on this method in Advanced Engineering Mathematics by Kreyszig
> in section Numerical methods for differential equations.
>
>
> I sent this uncleaned code to the Terragen guy and he made it a new
filter.
>
> ------------------------------------------------------------------
> LevelConnector uses these pseudo functions to communicate with the
> rest of the app.
>
>   SetHeight(x, y, hv);
>   hv=GetHeight(x, y);
>   Sel=GetSelValue(x, y); // is pixel selected
>   PixV=HVToFloat(hv);    // convert from internal type to float
>   hv=FloatToHV(PixV);    // convert to internal type
>
>
> file://def for the compare func used by bsearch
> typedef int (*fptr)(const void*, const void*);
>
> int compare(const int *a,const int *b)
> {
>   return(*a - *b);
> }
>
>
> // The main calculations. Processes data in u[].
> // Solves this Laplace's 2nd order differential
> // equation using Gauss-Seidel iteration.
> //
> //   d2u   d2u
> //   --- + --- = 0 + (force, not used)
> //   dx2   dy2
> //
> // Can be solved by FFT also but this iterative method lets user control
> // the accuracy of the solution. If all pixels are selected then this
works
> // like smooth with variable radius. Basically it just averages 4 adjacent
> // pixels into the middle pixel :) Note however that one iteration pass
uses
> // uses 2 values per pixel already computed in that pass so it's a kind of
> // cumulative smooth.
>
> void SolveEquations(int **s, float *h, float *u, int Nodes, int maxiter)
> {
>   int i,c,dir=1,Node;
>   float sum;
>
>   for(i=0;i<maxiter;i++)
>
{for(Node=(Nodes-1)*(1-dir)/2;dir*Node<=(Nodes-1)*(1+dir)/2;Node=Node+dir)
>         { sum=h[Node];
>           for(c=1;c<=s[Node][0];c++)
>             sum=sum+u[s[Node][c]];
>           u[Node]=0.25*sum;
>         }
>       dir=-dir;  // toggle direction
>     }
> }
>
>
> // Prepares data for SolveEquations
>
> void LevelConnector(int width,int height, int maxiter)
> {
>   int xo[4]={1,-1,0,0},yo[4]={0,0,1,-1};
>   int c,Nodes=0,*NodeP,Node,PixN;
>   int x,y,xt,yt,i,j,*Pos,**s;
>   float *h,*u,Sel,PixV,k;
>
>   // Count selected pixels
>   for(y=0;y<height;y++)
>     {for(x=0;x<width;x++)
>        {if(GetSelValue(x, y) == SELECTED)
>           Nodes++;
>        }
>     }
>
>   u=(float *)malloc(Nodes*sizeof(float)); // This is the actual height
data
>   h=(float *)malloc(Nodes*sizeof(float)); // aux data for equations
>   Pos=(int *)malloc(Nodes*sizeof(int));   // Stores pixel positions of
nodes
>   s=(int **)malloc(Nodes*sizeof(int *));  // aux data for equations
>
>
>   // clear arrays
>   for(Node=0;Node<Nodes;Node++)
>

>     s[Node]=(int *)malloc(5*sizeof(int));
>     if (s[Node] == NULL) {perror("memory error");exit(1);}
>     for(j=0;j<5;j++)
>       s[Node][j]=0;
>     u[Node]=0;
>     h[Node]=0;
>     Pos[Node]=0;
>   }
>
>
>   // Fill position table
>   Node=0;PixN=0;
>   for(y=0;y<height;y++)
>     {for(x=0;x<width;x++)
>        {Sel=GetSelValue(x, y);
>         if(Sel == SELECTED)
>           {Pos[Node]=PixN;
>            Node++;
>           }
>         PixN++;
>        }
>     }
>
>
>   // Make equations
>   Node=0;PixN=0;
>   for(y=0;y<height;y++)
>     {for(x=0;x<width;x++)
>        {if(GetSelValue(x, y) == SELECTED)
>           {k=0;c=0;
>    for(i=0;i<4;i++)
>              {yt=y+yo[i];
>               xt=x+xo[i];
>               if (yt == -1)      yt=height-1;
>       if (yt == height) yt=0;
>       if (xt == -1)      xt=width-1;
>       if (xt == width)   xt=0;
>               hv=GetHeight(xt, yt);
>               PixV=HVToFloat(hv);
>               if(GetSelValue(xt, yt) != SELECTED)
>                 k=k+PixV;
>               else
>                 {c++;
>                  PixN=yt*width+xt;
>          NodeP=(int *)bsearch(&PixN,Pos,Nodes,sizeof(int),(fptr)compare);
>                  if(NodeP == NULL)
>                    {perror("error in bserch");
>                     exit(1);
>                    }
>                  s[Node][c]=(int)(NodeP-Pos);
>                 }
>      }
>            hv=GetHeight(x, y);
>            u[Node]=HVToFloat(hv);
>    h[Node]=k;
>    // h[Node]=k-force;  past experiment
>    s[Node][0]=c;
>    Node++;
>   }
>        }
>     }
>
>   SolveEquations(s,h,u,Nodes,maxiter);
>   free(h);
>
>   // Write data in u[] back to hf
>   for(Node=0;Node<Nodes;Node++)
>

>       free(s[Node]);
>       hv=FloatToHV(u[Node]);
>       x=Pos[Node]%width;
>       y=Pos[Node]/width;
>       SetHeight(x, y, hv);
>     }
>
>   free(s);
>   free(Pos);
>   free(u);
> }
>
> K.K.


   Aaaaaaaargh!!!  You've frightened the s**t out of me now
  Kari!!  Whats a man to do?? (Need new specs after that! Me 'eads gone all
dizzy!! ;o)

                                ~Steve~

    PS, how do I save an image as .jpg, *in POV-Ray*, to 'send' to you very
'NICE' people to perooooose over? Simple question, I know, but I'M sure that
you can handle it, Kari............thanx.

Post a reply to this message